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Abstract 

In this paper we extend a method for iteratively improving slow manifolds so that it 
also can be used to approximate the fiber directions. In its original form the method was 
previously used succesfully by the first author and C. Wulff to obtain slow manifolds, 
including normally elliptic ones in Hamiltonian systems, with exponentially small error- 
field in analytic systems. The extended method is applied to general finite dimensional 
real analytic systems where we obtain exponential estimates of the tangent spaces to 
the fibers. The method is easily implemented numerically, which we demonstrate on 
the Michaelis-Menten-Henri model. Finally, we extend the method further so that it 
also approximates the curvature of the fibers. 
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1 Introduction 



Singularly perturbed systems involving different scales in time and/or space arise in a wide 
variety of scientific problems. Important examples include: meteorology and short-term 
weather forecasting [TTl UHl EH] , molecular physics and the Born-Oppenheimer approximation 
[2T] . chemical enzyme kinetics and the Michaelis-Menten mechanism [22], predator-prey and 
reaction-diffusion models [23], the evolution and stability of the solar system [151 [16] an d the 
modeling of tethered satellites [291 EO] - These systems can also be "artificially constructed" by 
a partial scaling of variables near a bifurcation [26] . The main advantage of identifying slow 
and fast variables is dimension reduction by which all the fast variables are "slaved" to the 
slow ones through the slow manifold. Dimension reduction is one of the main aims and tools 
for a dynamicist and the elimination of fast variables is very useful in for example numerical 
computations. In fact, since fast variables require more computational effort and evaluations, 
this reduction often bridges the gap between tractable and intractable computations. An 
example of this is the long time (Gyears) integration of the solar system, see [T5J [16]. See 
also [2] for a numerical treatment of slow-fast systems. 

Slow-fast systems and singular perturbation theory. We consider slow-fast sys- 
tems of the form 

d t x = eX(x,y), (1.1) 
d t y = Y(x,y), 

with a small parameter e. The vector-fields X and Y will in general also depend upon e, 
but we shall throughout suppress this dependency. In fact, one of the advantages of the 
method we use is that it does not require any smoothness of X and Y as a function of e. 
We assume that d y Y is invertible on the set {Y(x, y) = 0} for e = 0. This is precisely 
the meaning of y being fast. Here d z is used to denote the partial derivatives J-^, and we 
will continue to use this symbol regardless of what object is being differentiated. By the 
implicit function theorem the set {Y(x, y) = 0} can therefore locally be represented as a 
graph M = {y = r)o(x)}, and we introduce (x,y) = (xo,yo + r}o(xo)) to transform these 
equations into 

d t x = eX {x ,y ), (1.2) 
d t yo = Y (x , y Q ) = po(x ) + A (x )y + Rq(x , y ), 

with 

Po = -ed x T]X(x , 770), Aq = d y Y(x , Vo) ~ ed x rjd y X(x , r] ), 

Ro( x ,yo) — 0(Vo) an d X (xo,y ) = X(x ,r] + y ). The manifold M is not invariant since 
Y Q \ yo= Q = po, but it is close to being invariant as p — C(e) is small. 
In slow-fast systems one often connects ( 11. 21) with the system 

d T x = X {x ,yo), (1.3) 
ed T y = p (x ) + A (x )y + R (x , y ), 
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related to (11. 2\i by the scaling of time r = et. If we formally set e = in (11.2p and (11. 3p . 
then two limit systems are obtained. In ( II. 3p . the formal limit is singular leading to the 
algebraic equation: A (xo)yo + Rq{xq,Uq) = 0. Note that A and Ro are here evaluated at 
e = 0. Due to the singular nature of this limit the theory is also referred to as singular 
perturbation theory. The set of points M = {y = 0} satisfying these equations is called a 
slow manifold in the sense of MacKay [20], Definition 1, (or a critical manifold [11]) and the 
corresponding system: d r Xo = X (xq, 0) is called the slow subsystem^ On the other hand, in 
(II. 2p . the formal limit leads to the fast sub-, or frozen, system: d t y = A (x )y + -Rot^O; Ho) 
with xo now considered as a parameter. In this system, M is a set of equilibria and the 
linearization d T Syo = Ao(xo)Syo about these determine the classification of the slow manifold. 
In particular, if yo = is an elliptic or hyperbolic equilibrium, then the slow manifold Mq 
is said to be normally elliptic respectively hyperbolic at the point (xq,0). If M is normally 
elliptic and the vector-field is real-analytic then there exists an almost invariant slow manifold 
M e nearby. By "almost" it is understood that the error field, that is the normal component 
of the vector field restricted to the slow manifold, is of order 0(e~ c ^),C > 0. If the slow- 
fast system is Hamiltonian then M e can be made symplectic on which a (formally) reduced 
Hamiltonian system can be defined. These are basically the main results of [31J. These 
results are obtained using the method of straightening out which we now explain. 

The SO method. We now describe the method of straightening out as it is presented by 
MacKay in [20]. We will henceforth abbreviate this method by SO. It is this method which 
we in this manuscript seek to extend. The method is iterative, considering normal forms 
of the form (11.21) at each step of the iteration. To complete the first step of the iteration, 
consider the equation Y = 0, with Y as in (11.21) . This gives, by applying the implicit 
function theorem, a solution y = r/i(x ) close to 

Vl w -Aq Vo, (1-4) 

since Ro = O The graph M\ = {yo = r)i(xo)} will be an improved slow manifold. 
To show that this is indeed an improved slow manifold, one straightens out the new slow 
manifold by introducing y\ through yo = yi + Vi- Then the equations become 

d t xo = eXi(x ,yi), d t y x = Yi(x ,yi) = pi(x ) + A 1 (x )y 1 + Ri(x ,yi), 

with 

Pi = -^d x r} 1 Xo(xo,T}i) i (1.5) 

and so formally p\ = 0(e 2 ), since rji = 0(e) (11.41) . which is the measure of the error-field, 
an improvement from (9(e) to C(e 2 ). Here Y\ = —ed x r] Xo(xo, yi + 771) + Yo(xo, y\ + 771) and 
X\ = X (xo, yi + 771). We use d x rj rather than d XQ r) to avoid clutter. Continuing in this way, 
at each step solving Yi(xo,yi) = for = r] i+ i(xo) and then setting yi = y i+ x +r] i+ i(xo), we 
obtain an improved error at the end of each step which is a (9(e)-multiple of the previous 
error leading to formal (9(e n )-estimates. The method, though viewed slightly differently, is 
actually identical to the method suggested by Fraser and Roussel [6], [27]. In [12] this method 

1 Often slow manifolds are also required to be invariant, but this is too strong a requirement in the 
normally elliptic setting. This is the motivation for MacKay's definition. 
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is also referred to as the iterative method of Fraser and Roussel. This can be realized by 
introducing the partial sum 77" = Y^l^iVi an d expanding F n _i(xo, r] n (xo)) as 

Yn-xixo, r] n {x )) = -ed x rj n _iX n _ 2 (xo, rin-i + Vn) + Y n _ 2 (x , ?7 n ,_i + 7] n ) 

= -ed x r) n -iX n _ 3 (x Q , r\ n ^ 2 + Vn-i + rj n ) ~ ed x r] n ^ 2 X n ^ 3 (x , r\ n ^ 2 + Vn-i + rj n ) 

+ Y n _ 3 (x , T]n-2 + Vn-1 +Vn) = ■■■ 

= -ed x rj n - 1 Xo(x ,v n )+Yo(xo,V n )- 

The equation: 

-ed^Xoixorf) +Y {x ,ri n ) = 0, (1.6) 

defines the nth step of Fraser and Roussel's iterative method in which one solves for an 
improved slow manifold t] n , see [6l[27] and [12] where an asymptotic analysis of the method is 
given. The reference [12] does, however, not obtain exponential estimates. When the method 
is viewed within MacKay's setting we can also realize that we can in fact allow A to be an 
unbounded operator: It is only necessary to assume that A (x)~ l is bounded, making the 
approach potentially useful for partial differential equations. Note that p\ actually vanishes 
at a true equilibrium where X (xo,?7o) = 0, and the improved slow manifold M\ = {y± = 0} 
therefore includes all equilibria near Mq. This property is preserved when using the method 
iteratively. 

For Neishtadt's Hamiltonian example 

H = ^x 2 + iy 2 + v + eyf(u), 

with f{u) = Yln°=i e ~ n sin(rax) and symplectic form u = dx A dy + e~ x du A dv, one can by 
introducing ( = x + iy construct the optimal slow manifold and show that it cannot be 
improved beyond an exponential estimate. See [8]. The almost invariance is therefore the 
best one can aim for in a general setting for normally elliptic slow manifolds. Normally 
hyperbolic slow manifolds, on the other hand, persist. Nevertheless, it is not in general 
possible to show convergence of the iteration. This is due to the fact that the persistent 
manifolds, as with center manifolds, are not in general analytic p~|. 

Normally hyperbolic slow manifolds and their fibers. When M Q is normally 
hyperbolic then Fenichel's theory [U [5] applies and there exists a perturbed, invariant slow 
manifold M e for e 7^ nearby. Moreover, the stable and unstable manifolds persist. To 
explain the latter, consider at e = the fast fiber F Z ° = {(a?o> Uo)\ \\yo\\ < A} based at the 
point zq = (xo,0). If the real parts of the eigenvalues of A are all negative, then Mq for 
e = is asymptotically stable and all solutions on Fq° contract exponentially toward the 
base point zq provided A is sufficiently small. By Fenichel's theory the fast fibers Fq° perturb 
to F£° forming an invariant family F € = \J Zoe M e F*° along which solutions contract to M e . 
The invariance of this family is understood in the following sense 

$ f (F £ Z0 ) C Ffo (2o) , 

where $q is the time-t flow map of (11. 2p . The motion of any point on F e z therefore decomposes 
into a fast contracting component and a slow component governed by the motion of the base 
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point Zb = TTfZ of the fiber. In the physics literature a fiber is also sometimes called an 
isochron [3]. 

The fiber projection 7Tj is smooth, and so locally there exists a transformation (u,v) h-> 
(xq, yo), which is 0(e) close to the identity, mapping (jl.2p into the Fenichel normal form [TT] : 

d t u — eU(u), d t v — V(u,v)v. 

These are the ideal coordinates for the description of the system near the slow manifold, 
in particular useful when solving the problem numerically; the slow manifold coincides with 
the zero level set {v = 0} and the fibers of the form _pj n *' ^ have been straightened out to 
{(u, v)\u = u*, \\v\\ < A}. We will approach this ideal by first constructing a transformation 
(x, y) (-)■ (xo, yo) so that the x-equation, up to exponentially small error terms, becomes 
independent of y to linear order: 

d t x = e(A(x) + 0(y 2 )) + 0(e- c ^), 
d t y = A(x)y + 0(y 2 ) + 0(e- c ^). 

Then the tangent space to the fibers at (x*,0) will almost coincide with {(x, y)\x = x*}. 
Later we will also seek to remove the terms that are quadratic in y. Similar ideas have been 
developed in [31 [25] for center manifolds near non-hyperbolic equilibria. 

When M is of saddle type, with a stable manifold W S (M) of dimension n s j and an unstable 
manifold W U (M) of dimension n u j (rif = n^+n^j), then Fenichel's normal form takes a slightly 
different form: There exists a transformation (u,v,w) n- (xo,yo)j with dim{t>} = n s j and 
dim{w} = rij, which is (9(e) close to the identity, mapping ( 11. 21) into 

d t u = e(Uo(u) + Ui(u,v,w)vw), 

d t v = V(u, v, w)v, (1.7) 
d t w = W(u, v, w)w. 

Here Ui(u,v,w) : {v} x {w} — > M. ns is a bilinear function of v and w. The slow manifold is 
then given by {v = 0, w = 0} with stable manifold {w = 0} and unstable manifold {v = 0}. 
The transformation may only exist in a small neighborhood of the slow manifold so in general 
we need \\v\\ < A v and \\w\\ < A w . 

Reduction methods. There are several alternatives to the method of straightening out 
for approximating slow manifolds. We name a few others: The intrinsic low-dimensional 
manifold (ILDM) method of Maas and Pope [19], the zero- velocity principle (ZVP) [7] [32]. 
and the computational singular perturbation (CSP) method initially due to Lam and Goussis 
[T3l [T4] . but later extended by Zagaris and co-workers [33]. The ILDM method is based on 
the Jacobian of the vector-field and partitions this at each point into a fast and a slow 
component based on spectral gaps of the Jacobian. The ILDM approximation to the slow 
manifold is then defined as the locus of points where the vector-field lies entirely in the 
slow subspace. In general, this only gives an approximation that agrees up to 0(e) [T2]. 
Nevertheless, the method is still quite powerful as it can be used in systems where a small 
parameter may not be directly available. In the ZVP method an 0(e n )- approximation to 
the slow manifold is obtained as the locus of points where the (n + l)th time derivative of 
the fast variable vanishes. This method has been used in an equation- free setting in [7]. 



5 



A distinctive contribution of the CSP method is that it also approximates the fibers in 
the sense that it provides a set of basis vectors spanning the tangent spaces of the fibers 
at the slow manifolds up to 0(e n ). We will in this paper show that it is also possible to 
approximate these tangent spaces by adding an extra step to the method of straightening out. 
The transformations we use have direct interpretations, and the method can be formulated 
in a way that is similar to how Fraser and Roussel presented the method of straightening 
out, only involving evaluations of the initial vector-field and its Jacobian matrix. Finally, the 
method can be further extended so that it also approximates the curvature of the fibers. In 
principle higher order effects can also be accounted for, but this introduces a certain degree 
of complexity. In this paper we will therefore focus most of our effort on demonstrating the 
first part of the method which seeks to estimate the tangent spaces of the slow manifold. 
Once this approach has been established and demonstrated on the Michaelis-Menten-Henri 
model, we will consider removing the part of the slow vector-field which is quadratic in 
the fast variable, hence approximating the curvature of the fibers. One of the reasons for 
choosing the Michaelis-Menten-Henri model as our example is that all the calculations can 
be done explicitly. But moreover, it also allows for comparison with the results in [33] from 
the application of the CSP method. 

The SOF method. In this section we shortly describe our method for approximating 
the tangent spaces of the fibers. We will refer to this method as the SOF method - the 
extra F has been added to SO to indicate that the approximation of the fiber directions is 
build in as an extension of the original SO method. The method is based on normal form 
computations, see e.g. [9] section 3.3. We start from the real analytic slow-fast system: 

d t x = eX(x , y) = e(A(x ) + p (x )y + T(x , y)), 
d t y = Y(x , y) = p(x ) + A(x )y + R(x , y), 

with p describing the error-field on {y = 0} and R, T = 0(y 2 ). We assume that there are n s 
slow variables x G M™ s and rif fast variables y G M n/ . Following the 0(e~ 1 ) applications of 
the SO method we can take p = 0(e~ c/e ) ^T], and for the purpose of obtaining exponential 
estimates we can therefore ignore this term completely. The aim is to introduce a succession 
of transformations of the form Xi = Xi + \ + e<pi(x.i + i)y formally pushing the term in e -1 <9fXj + i 
which is linear in y to consecutive higher order in e. Let us consider the first step, introducing 
xo = x\ + €(po(xi)y so that 

d tXl = J' 1 (eA(x ) + e {e<9 x A0 o + - <p Q A} y + eO(y 2 )) 
= e (A + {ed x K<Po + IM> ~ M V ~ ed x (f> Ay + 0{y 2 )) (1.8) 

where J = I s +ed x <j)oy, I s = identity G ]R nsXri£! ; is the Jacobian of the transformation x\ i— > Xq, 
and where we have used the identity 

J' 1 = I s - ed x (f)oy + r l \ed x <i) yf '. 

The term in (11.81) which is linear in y is due to two contributions. The first one is due to the 
expansion of X(x ,y) — (pY(x ,y) in y, the curly bracket in (11. 8p . while the second one: 

Pi = -edxfaA, (1.9) 
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comes from the inverse of the Jacobian. Here d x (poA is understood column-wise: 

cy> A= (5,(0)^- --^(0)"/ A), 

(4>y = (4>y(x ) G M™ s being the ith column of 4> = 4>(x ) G M nsXn ^. We let O be the solution 
to the linear equation obtained by setting the first contribution, the curly bracket in (\1.8\i , 
to zero: 

ed x A(t> + no - <j) A = 0. 

This equation has a solution O close to /i A _1 , and the new error term fx\ P-9P ? which by 
construction is the only remaining term in e~ l d t X\ linear in y, is therefore formally smaller 
than the old error fiQ. There is an improvement from 0(1) to 0(e). Note also that 

Ai = A, A x = A. 

We will use these types of transformations successively in the proof, pushing the error term 
to higher order in e. One of the main results of the paper is that eventually the error is 
exponentially small: fx = 0(e~ c ^ e ). 

We present the first result formally in Theorem 12.11 which we prove in section [31 In 
section [4] we demonstrate how the method can be used computationally on the Michaelis- 
Menten-Henri model. In section [5] we present a result, Theorem I5.1[ on approximation 
of the curvature of the fibers. Theorem 15.11 excludes normally elliptic slow manifolds and 
neutral saddle-type slow manifolds where both A and —A, Re A ^ 0, are eigenvalues of A. 
This requirement appears in the construction of the appropriate transformations, where we 
encounter linear matrix equations of the form: 

A T f + fA = Q\ 

for the unknown matrices Solutions of this linear problem exist and are unique if and 
only if <j(A) fl a (—A) = 0, see [10], Theorem 4.4.6. Here o~(A) denotes the spectrum of the 
matrix A. The case where both A and —A, Re A ^ 0, are eigenvalues of the A leads to small 
divisors, as in the problem of analytic linearization [9], and it will be considered separately 
in section [6l 

As opposed to [31] the applications we have in mind are primarily for normally hyperbolic 
slow manifolds, where the fibers provide the directions of the stable and unstable manifolds 
along which the solutions relax to respectively escape the slow manifold. However, our results 
in Theorem 12. II still hold true for the normally elliptic case by providing coordinates in which 
the slow dynamics become almost independent of the fast variables to linear order. The 
improved slow manifold includes all nearby equilibria and, what is also new, the Jacobian of 
the vector-field at an equilibrium takes a very suitable form with the linearized slow dynamics 
being exactly independent of the fast variables. 

Notation and preliminaries. Let (X, || • and (y, \\ ■ \\y) be real Banach spaces, and 
Xc = X © iX respectively y<c = y © iy their complexifications with norms \\xi + ix2\\x c = 
11^111^ + H^Ha" and 1 1 2/1 + %2||y c = H2/2 ||y- Here we are primarily thinking of X = M. ns 

and y = M. nf , but other spaces will also be used and therefore prefer the general setting. 
As the proofs involve many transformations, we prefer to reserve the use of subscripts to 
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keep track of the steps in these coordinates changes, and will therefore use the alternative 
notation (x)i, 1 < % < n, for the ith component of a vector x G M. n . 

We will from now on denote all norms, including operator norms, by || • ||. We hope 
that this will not cause unnecessary confusion. Hopefully it will be clear from the context 
what norm is used. Then / : Uc — > 34: > with Uc an open subset of Xq, is analytic if it 
is continuously differentiable. That is if there exists a continuous derivative d x f : Uc — > 
L(Xc, yc), the Banach space of complex linear operators from Xc to equipped with the 
operator norm, satisfying the following condition 

\\f( x + h)-f(x)-d x f(x)(h)\\ = 0(\\h\\ 2 )- 

By real analytic we will mean analytic and real when the arguments are real. The higher 
order derivatives can be defined inductively and d™f becomes a map 

d:f:U c ^L n (X c ,y c ), 

from Uc into the Banach space L n (Xc, 3^c) of all bounded, n-linear maps from Xc x ■ • • x Xc 
(n times) into yc- See [23] Appendix A for a reference on analytic function theory in Banach 
spaces. 

When U is an open subset of X then we define U + i% to be the open complex x~ 
neighborhood of U: 

U + ix = {x E X c \d Xc (x,U) < X}, 

where dx c is the metric induced from the Banach norm || ■ ||. 
We frequently need the following Cauchy estimate: 

Lemma 1.1 \2J$ Assume that f : Uc — > yc is analytic and that f is bounded on the Xc-open 
ball B^(x ) C Uc- Then 

\\a ft mi ^ ^xeB^xo) \\f( x )\\ . . 

\\d x f{x )\\ < . (1.10) 



Remark 1.1 Consider / : U + %x — >■ 3^c analytic and bounded. Then we can apply this 
estimate to any xq G U + i(x — to obtain: 

sup ||a,/ M ||< sup ^ ll/W11 , 

which we will write compactly as 

ll/(s)llx 



This is the form of Cauchy's estimate that we will be using. Similarly, we will by || ■ \\ Xjl/ 
denote the sup-norm taking over the domain (U + ix) x (V + iv) of (x, y). 

Note also that the norm on the left hand side of ( 11.101) is the operator norm on L(Xq, yc) 
of complex bounded linear operators, while the norm on the right hand side is the norm on 

yc- □ 
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Remark 1.2 We write a m-linear form such as d™f(x) G L m (Ac,34c) evaluated diagonally 
on h e X as d™ f(x)h m . With this notation Taylor's formula reads: 

f(x + h) = f{x) + d x f(x)h + ■■■ + -J—dZ^fWh"- 1 

[n — 1)1 

r l (i _ s )n-l 
+ / { -^jy-d:f(x + sh)h n ds, 

with the integral remainder being bounded by — p sup 0<s<1 

2 Main result 

We consider the real analytic slow-fast system ( 11 .21) in the form 

<9tx = eX (x , y ) = e(AoOo) + fi (x )y + T (a;o, y )), (2-1) 
5*2/0 = *o(s , J/o) = Pot^o) + A)Oo)yo + #o(£o, yo), 

with n s slow variables and nj fast ones so that xq ElA + ixo C Ac = C ns and yo G V + zVo C 
y c = C Uf . Here U C X = R ns and V C y = R nf are real open subsets. Assume also that 
l|Po|| Xo <C(e). 

Theorem 2.1 Fix < x < Xo < < z^o- 27ien i/iere exists an e so i/iat /or a// 
e < e the SOF method constructs a transformation (x, y) (->■ (x , yo) which is e-close to the 
identity from {hi + i%) x (V + zV) to (U + ixo) x (V + iu Q ) mapping \2. 1\) into 

d t x = e(A(z) + fi(x)y + Q(x)y 2 + C{x, y)) (2.2) 
d t y = p(x) + A(x)y + R(x, y) 

with fi and p vanishing at true equilibria, 

l=\\p\\ K ,5 = \\p\\ K <0(e^), 

and T(x, y) = Q(x)y 2 + C(x, y), C = 0(y 3 ), 

||A — Ao||p ||^4 — ^ollx' \\T — ToW^v, \\R — Ro\\x,u < c 2^, 

for some constants c\ and c<i- □ 

We highlight that the estimates are not uniform in x an d \L- We also have the following 
corollary which provides a convenient form for the transformation in Theorem 12.11 

Corollary 2.1 If the eigenvalues of A Q all have non-zero real part, then there exists an eo 
so that for all e < e there exists a slow manifold M e of U.2\) and N e = 0{e~ x ) G N so that 
M e is given as the graph 

y = V (x ) + O(e- c ^), 
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with r] = J^n^iVn, where the partial sums rf = Y^i=iVi satisfy U.6\) for 1 < n < N e using 
the convention rf = 0. Furthermore, the tangent space of the fibers J 7 * at the base point 
z o = (xo,v(. x o)) i s given as 

^-r<( //+ ;^U )) +o(c " c,/,) )' (2 - 3) 

Here If = identity G M n / X ™/ and (f> = X^io^« cohere the partial sums <fi n = X^Lo'M^o) 
satisfy 

e(d x X + d y X d xV )r - ed x (j) n - l X Q + d y X - <j> n {-ed x rid v X Q + d y Y ) = 0, (2.4) 

for < n < N e using the convention <fi~ l = 0. The functions X , d x X , d y X Q and d y Y are 
all evaluated at (xq,t)(xo)), 

Proof Here Fenichel's theorem applies. The first part of the corollary then follows di- 
rectly from [31]. For the second part, note that each (p n solves (13 .4p below. Here A(xo) = 
Xo(xo,v(. x o)) an d A(xo) = —ed x r]d y X (xo, r]) + d y Yo(xo,r]) from which (12.41) follows by sum- 
mation over n. Also since the method generates a transformation of the form 

x = x + e<p(x)y + C(y 2 ), y = y + rj{x ). 

we obtain a tangent vector to the curve 6 = 6{{y)i) at (xo, v( x o)) as 

^ (0)= (e* + e ( ^K0) 4 )' 

Also (ei)j = 5ij Kronecker's delta, and (0) J = (0)*(xo) G W ns is the ith column of = 0(x) G 
l" sX "/. ■ 

We believe that these results, in particular in the form presented in Corollary 1, are useful in 
computations as the approximation of the relevant objects, the slow manifold and its tangent 
spaces, only require evaluations of the initial vector-field and its gradients. In particular, we 
believe that the approximations of the tangent spaces can be usefully applied in examples 
with many fast degrees of freedoms where one is faced with having to trade off accuracy 
with minimizing computational effort. From a given initial condition (xo,yo), near the slow 
manifold, one can approximate the fiber projection iif : (xo,y ) (£(,, ^(a^)), onto the base 
point, by solving the equations 

x = < pp + e<f>(xl™)y, (2.5) 

yo = y + v( x o), 

for y and x^ pp . The second equation gives y = yo—r](xo) which inserted into the first equation 
gives a non-linear equation for x^ pp . The right hand side of this equation is, however, e- 
close to the identity. In [33] it is stated that this projection is only 0(e), and therefore 
asymptotically in e not better than the "naive projection" (x ,y + r](x )) >->■ (x ,r)(x )). 
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However, this estimate is for fixed y. We believe it is more appropriate to highlight that the 
error is of the form: 

\\7r f (x ,y ) - (x^M^rm = O(ey'). 

In particular it is exact (up to the exponentially small error terms) if the tangent space is 
a hyperplane (which [33] also highlights). The different projections are illustrated in Fig. [TJ 
The linear projection accounts for the "initial slip" [3] along the slow manifold. 




Figure 1: Illustration of the different projections: naive, linear and exact. The linear projec- 
tion (x ,y ) i — v (%l PP , v( x b PP )) i s gi ven by the equations in (12. 5p . 

By approximating the fiber projection we can compute approximations to the dynamics 
having only to propagate initial conditions on 0(1) time scales, splitting the problem into 
first propagating the base point Xb = Xb(r), r = et, through 

d T x h = X (x,7](x b )), 

and then follow this by propagating y = yo(t) through 

d t yo = Y (x ,yo + r}(x )), 

using xo = Xb + e<fi(xb)y and the solution Xb = Xbij) obtained from the first step. We aim to 
investigate our results within this setting in future research. 



3 Proof of Theorem 2.1 



We first make use of the result from [31], to transform (12. ip into 

d t x = eX(x , y) = e(A(x ) + /J, (x )y + T(x , y)), 
d t y = Y(x , y) = p(x ) + A(x )y + R(x , y), 
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defined on the domain (xo,y) G (U + ix) x (V + iv) with \ — (x + Xo)/2 and ^ = u, and 
where p = 0(e~ Cl ^ e ). In fact, we initially ignore this term setting p = 0. Furthermore, X 
and Y are e-close to X respectively Y being given by 

X(x Q ,y) = X (x ,?7(xo) + y), (3.1) 
F(x , y) = -9 x ?7Xo(xo, 77(0:0) + 2/) + Fofao, 77(2:0) + y). 

Let if and C' A be so that ||A _1 || X < f and ||^.A|| X < C' A . 
We define the error 70 by 

7o = ||/A)|| X > 

and apply the transformation 

xq = xi + e0 o (a;i)y, 

where 0o solves 

e^A0 o + /io - <M = 0. (3.2) 

Cf. (jl.8p this transforms the system into 

= e(A(xi) + nxy + Ti(xi,y)) 
<% = + i?i(xi,y), 

with 



Remark 3.1 The modified function /ii vanishes at an actual equilibrium where A = 0. 
This implies that the linearization in these coordinates takes a very suitable form with the 
linearized slow dynamics dtSxi = ed x A5xi exactly independent of the fast variables. □ 

Note that 

\\xi - x \\ X: „ = e\\<t>oy\\ Xtl , < e7 cr, (3.3) 
where a = sup yeV+il/ \\y\\. From the linear equation ( 13. 21) we immediately obtain 
Lemma 3.1 If e < 1/(KC' A ) then the solution of A3.2fy satisfies 
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Proof Take 0° = jiA 1 , let r = y7 and introduce 0o = 0o + z so that (I3.2|) becomes 



2 

z = F(z) 



where F(z) = ed x A(<j)Q + z)A 1 . We have 

||F(z)|| x <eKC' A r<r, 

IM x <^<i. 

Here we have used the assumption e < 1/(KC' A ). The function F is therefore a contraction 
on B r C V + if and there exists a unique solution of f ]3.2f) with 

||0o|| x <2r = K 7o . . 

Using this lemma we can then estimate the new error using a Cauchy estimate 

7i = ll^i llxi ^ e — 7o, 

where Xi = X — £o- We now use this result successively, introducing x n = x n+ \ + e(fi n (x n+ i)y 
with (fi n solving 

ed x A(fi n +{J, n - <fi n A = 0, \i n = -ed x (fi n _iA, (3.4) 

on i 6 W + «Xn, Xn = X ~ YH=o Cnj for each n > 1. We take £ n = £ = 2eKC A at each step 
so that 

7n+l < ^7n < \ln < 2~ ( " +1) 7o, 
Sn ^ 

with 7 n = 1 1 /i n 1 1 xn 5 Xn — X ~ n £- Note also that 

n— 1 n— 1 

ll^n - ^ollxn,^ < ~ ^Hw ^ ecr 2 _i 7o < 2e<T7 . 

i=0 i=0 

Setting xtv e = X we realise that we can take N e = = = 0(e" 1 ) steps so that 

xo-x ' 



In. < 2 \ 4 ' kc aJ 1o . 

Including p = 0(e~ Cl ' e ) does not change the conclusion. The result therefore follows. 

Remark 3.2 Here we consider a fixed number applications of the SOF method, and show 
that the extension should only be iterated as many times as the first part has been iterated 
for the approximation of the slow manifold. To show this, we fix k G No, and apply the SO 
method k times to ( 11.21) so that the equations for (xq, y k ) = (xq, yo — V h ( x o)) are 

d t x = eX (x ,y k + r) k ), 

dtVk = pk(xo) + A k (x )y k + R k (x ,y k ), 
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with pk = 0(e k+l ). The SO method has then in fact been applied fc+1 times to the original 
equations (11. ip . Next, we apply the SOF method to these equations by introducing the 
transformation xq = x n+ i + ecj) n (x n+ i)y k where n = 0(1) solves (12. 4 p with r\ replaced by rj k : 



d t x n+ i = eX n+1 (x n+1 , y k ) = eyX - <p n p n 

+ [e(d x X + d y X Q d x r] k )r - ed x cj) n - X X Q + d y X - <f> n (-ed x r) k d y X + d y Y ) \y k 
- e(d x {r - <T x )Xo + <P n d xPk Ay k + 0(y 2 k ) 



= e [X - 4> n Pn + (ed x <f> n X + ^ y k + 0{y\ 

Here n = n — n ~ 1 = 0(e n ) and X , Y and their derivatives are all evaluated at (xo, r] k (xo))- 
Therefore the error, that is the term in X n+ \ linear in y k is formally of order 

ed x <p n X + erd xPk r = 0(e n+1 ) + 0{e k+1 ) = (e m ' m ^ +1 ), 

and, as expected, there is no improvement for n beyond k. A similar result holds true when 
considering the transformations in section [5] that seek to remove the terms in the slow vector 
field that are quadratic in the fast variables. □ 



4 Michaelis-Menten-Henri model 

In this section we demonstrate our method on the Michaelis-Menten-Henri model 

d t x = eX (x , y) = e(-x + (x + k - X)y), 
d t y = Y(x , y) = x - (x + 

for enzyme kinetics [31] . Here Xq and y are non- negative concentrations and the parameters 
satisfy k > A > and < e < 1. Setting Y(x ,y) = gives y = r] (x ) = an( j so 
U — Ho + ^oO^o) transforms the system into 

d t x = eX (x , y ) = e ( + (x + k - X)y ) , (4.1) 

V. x + k ) 

dm = Y (x , y ) = ^-^e - ^ + k + {xq + k)2 ej y , 
Therefore if k ^> e, so that 

A = d y Y (x Q , 0) = x + k + £K ( X ° + h A ) £j 

[Xo + K J 

Xo being non-negative, then the system is slow- fast with n s — 1 and rif = 1. The variable Xo 
is slow with <9 4 x = 0(e) and y is fast with (c^Yq) 1 = 0(1). 
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Analytic expressions of 77 and <j) to 0(e 2 ) 

We now obtain analytic expressions for 77 and <fi. Since the model is actually linear in the fast 
variable the SOF method only involves the solution of linear equations. First, we introduce 
771 satisfying Y (x ,r]i) = 0: 

(x + k)((x + k) 3 + €k(x + hi — A)) 

Then we define Yi(xo,yi) = —d x r] Xo(xo,r]i + y\) + Yq{xq,t}\ + and determine 772 from 
the condition Y\ (3:0,772) = 0. We obtain 

^2 = ~ 7 e + O e , 

(x + «) 

and therefore 

2/o = »7 =Vi + V2 = - 4 e ; 7 ^7 e + ( e ) > 

(x + «) (x + K ) 

as a 0(e 2 )-approximation of the slow manifold.- The error-field is 

/ \ . X 3 kx q (k 2 — 12kxq + 15 Xq 2 ) o m , 4X . 
p(x„) = Y 2 (x , 0) = ± ^e 3 + O (e 4 ) . (4.2) 

(X + K) 

To approximate the fiber directions we introduce y through 7/0 = rj 2 + y and compute 

Axq (k — A + Xq) k Xxq 



A = X (x ,v (xo)) = , 

(k — A + Xq) k\x (k 2 — 2kA + (k + 3 X)Xq) 
{k + x )' 



-€ 



2 



(e 3 ) (4.3) 



A = W zo, . 2 (x )) = -«-«,- ilipAl^l e _ («- 3s ) (*- A + *,,)* 

+ O (e 3 ) , (4.4) 

HO = X + K — \. 

Inserting this into (13. 4p with n = we obtain 

fl Xq + K — A (x + K — A) K (Xq — 2 A + k) rn( 2> 



A - e<9 z A x + « (x + k) 4 

At the next step we first compute the new error 

AH = -ed x ^A = "-^V + ( e ') 
(« + ^o) 

and via (13.4p with n = 1 we solve for 0x 

A*i A 2 z 



-e + 0(e 2 ). 



A - ed x A ( K + x ) 



-e + 0{e 2 ). 
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Then 



x (ft - 3x ) A 2 , 3 . 

^2 : — ; ^ e + <-A e J> 

(« + ^0) 



so that y i a fl3.4[) with n — 1 becomes: 



x (k — 3 £0) A 3 
Oo + «0 ? 



e 2 + G(e s 



Let 



+ 0i + fa- 



, 2 x + ft — A (ft 3 — 3ft 2 A + 2x ft 2 — 3 k Ax + x 2 n + 2 k A 2 + X 2 x ) . ,_ 7 

4 = ; + ; zz e - (so + «) 

^0 + « (x + ft) 

ft 2 (-A + k) (6 A 2 - 6 k A + ft 2 ) + ft (-A + ft) (-2 A + ft) (-2 A + 3 ft) x 



+ (-3 A 3 - ft 2 A + 3 ft 3 ) x 2 + ft (ft + 3 A) x 3 ^ e 2 + C(e 3 ). 



Then the span of the vector 



1 + ed x r] 




1 



(XQ + K-A) 
(XO + K — A)k 



+ 



/ (re(-A+ft)(-2 X+k) + (2k-\)(~X+k)x +kx 2 ^ \ 
(xo+re) 

(re (-A+k)(k--3 A) + (2re 2 -re A-2 A 2 )x +(k+3 A)x 2 )k 



V 



G(e s 



(x +re) 



gives a (9(e 3 )-approximation of the tangent space. We have left out the complicated 0(e 3 )- 
terms. The transformation (x,y) 1— > (xoj2/o) = (x + e<p 2 y,y + r] 2 (x + e</> 2 ?/)) transforms the 
system into: 



_ (x + ft ^ e + 0(e 2 ) ) j/2 + 0(e 2^ 



(x + ft)' 

% = p(x) + + 0(e 4 )) y + 0(q/ 2 ), 



(4.5) 



with 



/'3 



xA 4 (ft 2 — 12ftx + 15 x 2 ) 3 4 
- ■ -g e + C(e 

[x + ft) 



and where p, A and A are given in ( I4.2p . ( 14. 3 p respectively (14. 4p . 



Numerical computations of 77 and 

In Fig. |2] (a) and (b) we have compared the solution (x &, yob) of (14. ip initiated at the base 
point (xq 6 ,^(xq 6 )) with (i) the solution (xq,?/q) initiated at (xq 6 ,^(xq 6 )) + (dashed) 
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and with (ii) (by the naive projection) the solution (xq,i/q) initiated at (xq 6 , v( x ob) + s ) (^ u ^ 
line) for different values of s and for a time integration of length t = e -1 . The parameter 
s measures the distance from the slow manifold and v is the tangent vector to the fiber 
at the base point (xq 6 , r)(xQ b )) determined through and the equation ( 12.30 . We have set 
k = 2, A = 1 and x® = 1.5, and have computed rj and </> numerically using the procedure 
described in Corollary 12. II We have used 4 iterations on both r\ and resulting in error-fields 
of ~ 10 -7 respectively ~ 10 -10 for e = 0.1. The comparison is made through 

Vi = ||(x 06 ,y &)(e -1 ) - (4,2/o)(e -1 )ll> 
and 

vu = IKxot^ofcXe" 1 ) - (x^y^ie-^l 

In (a) e = 0.1 while e = 0.01 in (b). We see that Vi <C va and compute Vi ~ C(s 2,009 ), 
whereas vu w (^(s 1 - 000 ). 




(a) e = 0.1 (b) e = 0.01 



Figure 2: The errors from approximating fibers using the tangent spaces (vf. dashed lines) 
and naive projections (vu'. full lines) as functions of the distance from the slow manifold. 
Here k — 2, X — 1 and the initial condition on xq is x^ b = 1.5. 

The errors in if 1 and </> n : 

£ v = sup | - d x r}X + Y \, 

X 

respectively 

£ = sup \e(d x X + d y X d x 7])(f) - ed x (j)X + d y X - (f)(-ed x r]d y X + d v Y )\, 

X 

for e = 0.1 are shown in Fig. [3] as a function of the iteration number. These are relevant errors 
since if S v = then y = rj is an exact slow manifold and if = then the transformation 
Xq = X\ + e<j)(xi)y removes the term in d t Xi that is linear in y exactly. The error in <j) and rj 
are observed to be the order of machine precision ~ 10~ 14 after 8 respectively 10 iterations. 
There is no or little improvement beyond this number. Note that this in agreement with the 
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estimate of N e = 0(e x ) ~ 10 for e = 0.1. It should also be mentioned that to approximate 
derivatives we use the five-point stencil: 

f{x) « ^{~f{x + 2h) + 8/(x + h)- 8f(x -h) + f{x - 2h)), 

the error being fg/^Oco) — Oi^), xq E [x — 2h,x + 2h). We have used h ~ 10~ 2 which 
gives an error of ~ 10 -8 . 

10- 2 

10~ 4 

-©- i<r 8 

■! ,o-° 

o 

1— 

LU ,„-« 
10~ u 



" 2 4 6 8 10 12 

Iteration number 

Figure 3: The errors E n (o) and (x) for e = 0.1 as a function of the iteration number. The 
functions rj and <fi are computed numerically using the recipe in Corollary 12.11 Here k — 2, 
A = 1 and the initial condition on xq is x^ b = 1.5. In the computations leading to Fig. [2] we 
have used 4 iterations on both rj and (p. 

In Fig. H]we have taken s = 0.5 and consider Vi and Va as functions of e. Again, we see 
that Vi <^ Vu and compute t>j ~ 0{e 2 ) whereas vu ~ C(e). This discrepancy is, however, 
exceptional as it is due to the fact that the Michaelis-Menten-Henri system is linear in the 
fast variable yo. 

Finally, solutions (xob,yofc) and (x l ,y l ) of (14.11) are shown in Fig. [51 The solution (x l ,y l Q ) 
is initiated on the fiber of the base point with x^ b = 1.5, at a distance s ~ 0.52 from the 
base point. The full line near y = is the slow manifold. ixob,yob) and (x l ,y l ) at 9 
different times are indicated by x respectively o's. For illustrative purposes we have chosen 
the relatively large value of e = 0.4. It is observed that, at least approximately, the solution 
(x l Q, t/q) contracts along the fiber directions indicated by the dashed lines moving from upper 
left to lower right. The fiber directions are approximated as hyperp lanes through <ft. 

5 Approximating the curvature of the fibers 

In this section we show how the SOF method may be further extended to also approximate 
the curvature of the fibers. According to Theorem 12.11 and Corollary 1 the transformation 

y^y = y + 7](x ), x i-> x = x + e(p(x )y, 
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Figure 4: The errors from approximating fiber directions using the tangent spaces (dashed 
lines) and naive projections (full lines) as functions e. Here s = 0.5, k = 2, A = 1 and the 
initial condition on x is = 1.5. We have used 4 iterations in the computations of rj and 
0. 




Figure 5: Solutions (#o&,2/o&) anc ^ ( x o>2/o) °f (14.11) for e = 0.4. The solution (x ,y l Q ) is initiated 
on the fiber corresponding to the base point with x^ b = 1.5, at a distance s ~ 0.52 from the 
base, and is observed to contract to the solution (xo&,2/o&) a l° n g the fiber directions. The 
fiber directions are indicated by the dashed lines running from upper left to lower right. 



transforms (JET]) into f l2T2|) : 



d t x 

d t y 



e(A(x ) + Qo(x )y 2 
A(x )y + R(x Q ,y), 



C{x ,y)), 



(5.1) 



with C = 0(y 3 ) and R = 0(y 2 ) up to exponentially small error. We can easily obtain an 
explicit expression for the quadratic term Qoy 2 , a vector of symmetric bilinear forms, in 
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terms of the known functions: X , Y , rj and 0, through the equation Xo = Xo + ecj)(xo)y. We 
will write the ith component of Qoy 2 as 

(Qoy 2 )i = {y, Q l y), i<i<n s , 

where Q l Q = Q l {x) is a symmetric nj x riy-matrix. Recall that we use the notation (z)i to 
denote the ith component of a vector z. Here we have also introduced the real inner product 
(a, b) = XTii( a M & )i- B Y introducing 

x = Xi + e^ (xi)y 2 , 

with ipo a vector of symmetric bilinear forms, we therefore obtain 

dtxi = (I s - ed x ^y 2 + J-\ed x ^ y 2 ) 2 )e(A + {ed^y 2 + Q y 2 - 2^o(y)(Ay)} + 0(y 3 )), 

(5.2) 

where J = I s + ed x ipoy 2 is the Jacobian of the transformation £i i— > x . Here if)Q(y)(Ay) is 
understood as 

The new error, that is the term in e~ l d t X\ which is quadratic in y, can again be decomposed 
into two separate contributions. One term comes from the expansion of dtXo — (d y (ipoy 2 ))dty, 
the curly bracket in (I5.2p . while the other one is due to the inverse of the Jacobian. As for 
the linear case, we choose the unknown function ipo so that the curly bracket in ( 15.2ft vanishes 
for all y. This gives 

e d {x)j (A)i^ + Q l o - A T 4 - ip l A = 0. (5.3) 
i=i 

By Theorem 4.4.6 in [10] this system has a unique solution ip l Q for e = iff cr(A) Dcr(— A) = 0. 
Therefore we must exclude the elliptic case and the neutral saddle scenario where both A 
and —A, Re A ^ 0, are eigenvalues of A. We consider this separately in section [6] below. Note 
moreover that by taking transposes: 

Qh - A T (%) T - m T A = 0, 

and so this solution is symmetric. The solution perturbs to a symmetric solution for e ^ 
but small; e Y^j=i ®(x)j (A)^ ^ s a ^ so symmetric. The solution satisfies 

||Vo|lx = max||^|| x <iqQo|| x , 

for some constant K depending on A" 1 . Then the new error becomes 

Qi = -ed x il> A, \\Qi\\ x -(< — - — HQollx) 

which also vanishes at exact equilibria where A = 0. As for the linear case we have that 
Ai = A and A\ = A. Using such transformations successively it is therefore possible to 
approximate the curvature of the fibers up to exponentially small error. Formally we proceed 
as in the proof of Theorem 12 . 1 ( the most important ingredient in the proof being the continued 
reduction of the domain together with the applications of Cauchy estimates to control the 
derivatives. 
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Theorem 5.1 Assume that the assumptions of Theorem \2. 1\ hold true and that we have used 
the SOF method to transform (OOP into 112.2}) . Assume furthermore that o~(A) ncr(— A) = 0. 
Fix < x < X- Then there exists an e"o < eo, where eo /rom Theorem \2.1\ and an N € = 

0{e~ 1 ) e N so that for all e < eo i/ie sequence of transformations x n = x n+ i + eip n (x n+ i)y 2 , 
< n < N e - 1, where ^ E M n / Xn / so/ves 

e d (x) . (AUi + Q l n - A T ^ n -tiA = 0, l<t< n„ (5.4) 

the quantity 

n 2 _ / given by Eq. ( 15. ip for n = 0, 
^ = I -e^^Ay 2 /orn>l, 

fremg i/ie term in the expression for dtx n which is quadratic in y, eventually transforms $2.2}) 
into 

dfx^ = C (A&0 + C(x^y)) + 0(e- 5 ^), (5.5) 
% = A(x^)y + R(x# t ,y) + 0( e -^). 

Here {x^y) E (U + z%) x (V + zV), C = 0(y 3 ) and 

II C — C|Ix,h> ll-R — R\\x,ll — ^e, 

/or some constants c\ and c^- Also the 0(e~ Cl ^ e ) error terms in / 1 5. 5}) vanish at true equilibria. 

The transformation (x,y) >->■ (x,y) — (x + eipy 2 ,y) with ip = Xli^o" 1 ^ differs from the 
product (xft ,y) >->■ ■••{x\,y = y) •->■ (£o>2/) 0(y 3 )-terms and the equations for (x,y) 
therefore takes a similar form to l}5.5}) : The set {y = 0} is almost invariant and the y-space 
provides an almost 0(y 2 )- approximation to the fibers. In terms of the (xo,yo) -variables this 
quadratic approximation, parametrized by y, takes the following form: 

xo = x + ecpy + eipy 2 -, 

Vo = y + V + d x r}(e(f>y + eipy 2 ) + ^<9 2 ?7(e#) 2 , 
for the base point (x,r](x)). □ 



Analytic expression of ip to 0(e 2 ) for the Michaelis-Menten-Henri 
model 

We now apply this principle to the Michaelis-Menten-Henri model. We start from (14. 5 p where 
the quadratic term in e~ 1 d t x is of order 0(e): 



Qi = J* + «-»\ + 0(e .). 

[X + K) 
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We have therefore denoted this term by Qi rather than Qq. Then ipi solves (15. 4p with n = 1 
and ipo = 0: 



Qi (x + k - A) A 2 



so that 



Finally 



_ A 2 (2(x + At)-3A)x 

2(X + K) 



X 2 (2(x + .)-3X)x 

A(x + k) 6 K 



Let %p 2 = + ^2: 

, 2 (x + At — A) A 1 



V> a = — — — + 7 (z + k)" 6 ( 2 « (-A + «) (8 A 2 - 9 k A + 2 At 2 ) 
2 (x + At) 4 V 

+ (38 k A 2 - 44 At 2 A + 12 At 3 - 7 A 3 ) x + (4 A 2 + 12 At 2 - 22 « A) x 2 + 4acx 3 ^ + C(e 3 ), 

then, in terms of the original (xo, y)- variables, we have obtained the following quadratic 
0(e 2 )- approximation of the fiber with base point (x,rj(x)): 



x = x + e(j) 2 (x)y + eip 2 {x)y 2 

( (Sc + k-X) k (-A + k) (k - 2 A) + (2 k - A) (-A + k) x + kx 2 9 

= 1 + ("StT^ + a liT? 6 + 0(6 ' 1 S 

V 2 (x + At) / 
2/o = y + r/ 2 + d x rj 2 {e(j) 2 y + e^y 2 ) + ^«9 2 7? 2 (e0 2 y) 2 



AtxA At A X (At 2 + XAt — 2 At A + 3 A x) 2 , 3. / (x + At — A) (At — 3x) At A 

[X + At) (X + At) V (X + At) 



c 2 



(x + At — A) (2 At — 3x) At A , 



parametrized by ?/. We have here chosen only to show the (9(e 2 )-terms. In particular, 
(x,y) = (x + etp 2 (x)y 2 , y) transforms (14.51) into: 

d t x = e (A(x) + to{x)y + Q 3 y 2 + 0(e 2 y 3 )) , 
d t y = p(x) + (A(x) + G(e 4 )) j/ + 0(ey 2 ), 

with 

V 2 (5 + At) 5 1 ^ 
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Numerical computation of ip 



Fig. [6] shows the error 

v = ll(z&,2/&)(e _1 ) - (S:,y)(e~ 1 )\\ 

with (x, y) being the solution initiated at points along the quadratic approximation of the 
fibers with base (xb,yb) obtained from numerically computing the ip^s, as a function of 
the distance s G [0.5, 10] (s ~ y) from the slow manifold. We have used 4 iterations in 
computing ip giving rise to an error of ~ 10~ 8 for e = 0.1. The error decreases as ~ 0(s 3m3 ) 
in agreement with the analysis. 



10" 




10 - 



Figure 6: The error from approximating the fibers in the Michaelis-Menten-Henri model 
using a linear (full line) and quadratic (dashed line) approximation of the fiber as a function 
of the distance s G [0.5, 10] from the slow manifold. Here k — 2, A = 1, e = 0.1 and the 
initial condition on x is x® = 1.5. The error from the quadratic approximation decreases as 
pa 0(s 3 ). We have used 4 iterations in the computations of rj, <p and i/j. 



6 Saddle type slow manifolds 

Theorem 15 . 1 1 excludes the case of a neutral saddle type slow manifold. We proceed, as always 
in normal form calculations, by solving for what can be solved for, and therefore relax the 
requirements for the transformation. As before we are interested in a procedure described 
by simple equations and therefore choose not to use Fredholm's alternative to split the linear 
homological equations. We only aim to provide an outline for such a procedure here. 
Assume that Theorem 12.11 has been applied and consider the equations 

d t x = e(A(x ) + O(f))), 
d t y = A(x )y + O(f), 
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with (xo, y) 6 {U + ix) x (V + iu) related to the old variables through the equations 

x = x + &f>(x )y, 
Vo = v(%o) + y, 

and where exponentially small terms have been ignored. We assume that A has eigenvalues 
with positive and negative real parts. Before seeking to remove terms in the equation for xq 
that are quadratic in y, we first seek a splitting of y into "stable and unstable parts". To do 
so we first perform a change of basis 



y = V(x )y , A = V~ X AV 



A s0 
A u0 



where a(A s0 ) C (— oo, — A], a(A u0 ) C [A, oo) for some A > 0, so that 
d t x = e(A + O(y 2 )), 

dm = {A + eB ) y + O(f ), B = -V^d x VA. 
We then apply transformations of the form 

where (pi(xo) G ~R nfXnf for xq G U, to remove the off block diagonal terms of the resulting 
BiS. Let us consider the first step. Introducing y = (If + e(p )yi gives 

dm = (I f + e Vo )-\d t y - e 2 d x <p Ay 1 + 0(yD) 

= {(// + e^o)" 1 (Aq + eB ) (I f + e^ )} ft (6.1) 
- e 2 (I f + etpo^dxtpoAy! + 0(y 2 ). 

Set Jf = If + e<^o and introduce the following splitting of -Bo 

d ( Bso B su0 \ 

B = ( ^ s0 ® 1 B = ( ^ ^ su0 
\ 5 u o/ ' \-B us0 

We will collect the block diagonal parts -B i and A into Ai = ( ^f 1 ^ ) = A + e£>oi 



, A ul/ 

with A s i = A s0 + eB s0 , A uX = A u0 + eB u0 , and seek to remove B 02 . For this we expand the 
matrix shown with curly brackets in (16.11) : 

Jj 1 (A + eB ) J f = (If - eif + e 2 J/Vo) (^0 + e^o) (// + e(p ) 

= A 1 + e {A^o - PoA) + 5 02 } (6.2) 
+ e 2 (B cp - <p B - cp (A + eB )ip + Jj l yl(A + eB )J f ). 
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We then pick tpo of the form 

= 

so that the curly bracket in ( 16. 2 I) becomes 



(fisuQ 
</?usO 



A„a \ fA s0 



Aso^suo — V^swo^itO + B su q 

AuO^usO — PusoAsQ + B us q 

which we set to zero: 

Aso^suo — <PsuoA u o + B su0 = 0, A u0 ip us0 — ip us0 A s0 + -B us o — 0. 

These equations are solvable for (p su0 respectively ip us0 since by assumption we have a(A s0 ) n 
ct(74 u0 ) = 0. Also the solutions satisfy 

<^l|5o 2 || x , (6.3) 

for some constant K. The new error 

.61=6 (B ip - ip B - ip (A + eB )ip + J~Vo(A) + tB ) J) - eJ~ r d x tp Q h, 

is 0(e). Through (16.31) and a Cauchy estimate for the derivate d x <po, this error can in fact 
directly be estimated by £ -1 C7||.Bo2||x 011 the smaller domain IA + i(x — f° r some constant 
C > 0. We then again argue that for e sufficiently small we will after N e = 0{e~ 1 ) steps of 
this procedure have that B^ = C(e _C//e ). 

We then continue by considering = (y s , y u ) and the equations 

d t x = e(A(x ) + Q s0 (x )y 2 + Q u o{%o)yl + Qus(xo)(y s )(yu) + 0(3), 

d t y s = A s (x )y s + O(2), (6.4) 

dtVu = A u (x )y u + 0(2), 

where 0(2) and 0(3) denote terms that are quadratic respectively cubic in y^ and where 
we have ignored the exponentially small terms. The eigenvalues of A s all have negative real 
parts while the eigenvalues of A u all have positive real parts. In accordance with (11.71) . we 
will now aim to remove the terms in the equation for xq that are quadratic in y s and y u \ to 
avoid small divisors completely we will not seek to remove terms of the form Q su {ys){yu)- 
We therefore apply the following transformation to the slow variables: 

x Q = xi + e^svyl + eipuoyl, (6.5) 

to push the current errors Q s q and Q u q to higher order in e. We obtain: 

dtxi = J _1 ef A + led x k(4i s0 y 2 s + 4> u0 yl) + Q s0 y 2 s + QuoPl 

- 2ip s0 (y s )(A s y s ) - 2ip u0 (y u )(A u y u ) \ + Q us (y s )(yu) + 0(3) J. 
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Here we have suppressed the dependency on x\ and by J = I s + ed x ip s oy1 + ed x ip u oyl denoted 
the Jacobian of the transformation X\ i— >■ xq. Again we set the curly brackets to zero: 

ed x A(ei; s oy 2 s + e^ u0 yl) + Q s0 y 2 s + Quovl ~ ^so(y s )(A s y s ) - iJ u o(y u )(A u y u ) = 

for all y = (y s ,y u ). As in (15. 4p . this gives 

e ^ d {x) . (A)^o + Qlo ~ A^lo - Tpi A s = 0, 1 < z < n s , 
j'=i 

6 ( A ) A + Qlo - - tioAu = 0, 1 < % < n s . 

These equations are solvable for e sufficiently small as by assumption cr(A s ) flcr(— A s ) = = 
cr(yl u ) D cr(—A u ). The new errors are therefore 

Qsi = -ed x ip s0 A, Q ul = -ed x ip u0 A, 

which are due to the inverse of the Jacobian. They again vanish at true equilibria where 
A = 0. We iterate the procedure and apply Cauchy estimates to conclude that the errors are 
eventually exponentially small. The stable fiber of the invariant slow manifold {y s = 0,y u = 
0} then coincides with {y u = 0} while the unstable fiber coincides with {y s = 0} up to terms 
cubic in y s and y u and exponentially small terms. These spaces can through the functions 
(pi, <fi, ip si and ip ui be described in terms of the original variables. 

7 Conclusion 

In this paper we developed a new method, the SOF method, as an extension of the method 
of straightening out (SO method) so that it can also be used to approximate fiber directions. 
The method is based on normal form computations. After having approximated the slow 
manifold using the SO method, the extended method constructs a transformation of the slow 
variables as a product of a finite sequence of transformations, each obtained as the solution of 
a linear equation, so that the slow dynamics becomes almost independent of the fast variable 
to linear order. See Theorem 12.11 This gives an approximation, with exponentially small 
error, of the tangent spaces of the fibers. The method is easily implemented numerically as 
it only involves evaluations of the initial vector-field and its Jacobian matrix. See Corollary 
1. The method was also extended further in Theorem 15.11 so that it approximates, again 
exponentially well, the curvature of the fibers. This result holds true even when the slow 
manifold is of saddle type. In the saddle case, we only require that it is not neutral in the 
sense that there does not exist a contraction and an expansion rate of equal magnitude. 
For this scenario we followed the general philosophy of normal form theory and solved for 
what could be solved for. In agreement with Fenichel's normal form, we were able to con- 
struct a procedure, valid for any saddle type slow manifold not just the neutral ones, that 
removed terms in the slow vector field that were quadratic in the fast variables associated 
with the contraction respectively the fast variables associated with the expansion from the 
slow manifold up to exponential small error. 
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